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ABSTRACT 


The current study extends the application of computational fluid 
dynamics to three-dimensional high-lift systems. Structured, overset grids are 
used in conjunction with an incompressible Navier-Stokes flow solver to 
investigate flow over a two-element high-lift configuration. The computations 
were run in a fully turbulent mode using the one-equation Baldwin-Barth 
turbulence model. The geometry consisted of an unswept wing which spanned a 
wind tunnel test section. Flows over full and half-span Fowler flap 
configurations were computed. Grid resolution issues were investigated in two- 
dimensional studies of the flapped airfoil. Results of the full-span flap wing 
agreed well with experimental data and verified the method. Flow over the wing 
with the half-span was computed to investigate the details of the flow at the free 
edge of the flap. The results illustrated changes in flow streamlines, separation 
locations, and surface pressures due to the vortex shed from the flap edge. 
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CHAPTER 1 
INTRODUCTION 

1.1 MOTIVATION FOR RESEARCH 

High-lift systems continue to be an important consideration in the design of 
transport aircraft. This is due to the tremendous impact they have on the 
performance of the aircraft. For example, when considering a large, generic, twin- 
engine transport a 1.5% increase in CLmax will increase available payload by 6600 lb. 
at a fixed approach speed [1]. Similarly, a 1% increase in take-off L/D will increase 
the payload by 2800 lb [1]. According to Butter [2], a 5% increase in landing CLmax 
is approximately worth a 25% increase in payload. These relationships illustrate 
how even a small increase in high-lift system performance can make a big difference 
in the profitability of an aircraft. 

High-lift system performance is not the only issue, however. Mechanical 
complexity also factors into the overall evaluation of the system. The trend over the 
years has been to resort to added complexity to gain high-lift system performance 
(Figure 1 [3]). This tends to reduce the increased payload benefits illustrated 
previously because of the added weight and maintenance costs of the complex high- 
lift system. It is possible to gain performance without the added complexity, as 
illustrated by the Airbus aircraft in Figure 2 [4]. The Airbus aircraft (A300-A340) all 
show a similar CLmax 1° the Boeing aircraft (B727-747) at a lower complexity level. 

A better understanding of high-lift flows is needed to continue increasing the 
performance of high-lift systems while decreasing their complexity. It is likely that 
the ability to analyze three-dimensional flows will play a key role in such an 
improved understanding. Butter predicts that improvements in maximum lift will 
come from the ability to design for high Reynolds numbers and from a better 
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understanding of three-dimensional flows [2]. This study addresses three- 
dimensional flows through the computational analysis of a simple high-lift system. 

1.2 PREVIOUS WORK 

The high potential for payoff associated with high-lift system aerodynamics 
has stimulated much research in this field. Almost all of the literature falls within 
two categories: survey papers or two-dimensional studies. This section of the thesis 
is dedicated to summarizing the major high-lift research to date. 

The standard for survey papers in the high-lift arena was established by A. 
M. O. Smith's classic Wright Brothers Lecture [5]. Smith's paper started the 
"modem" era of high-lift system understanding. He dispelled the popular belief 
that the elements of a high-lift airfoil worked together by merely providing 
boundary layer control via the slots between elements. He was the first to illustrate 
the potential flow relationship between the elements of a high-lift airfoil. Since then, 
there have been a number of other survey type papers on the subject [1-4, 6-8]. All 
of these build on the findings of Smith [5] to identify areas of high-lift flows that are 
not fully understood. As an example, Meredith [1] summarizes some "viscous 
phenomena affecting high-lift system performance" as: 

• Attachment line transition from laminar to turbulent 

• Relaminarization of turbulent boundary layers 

• Transition of boundary layers 

• Viscous wake interaction 

• Confluent wakes and boundary layers 

• Shock/boundary layer interaction 

• Separated flows 



These survey papers also typically discuss the current state-of-the-art analysis 
techniques, generally commenting on the lack of a suitable analysis tool. A final 
element that is common throughout most survey papers is a summaiy of past high- 
lift research at major aircraft companies such as Boeing [6], Fokker [8], or British 
Aerospace [2]. In all, these papers offer a cursory description of the flow physics 
associated with high-lift systems, in addition to identifying problem areas and 
possible analysis methods. 

The vast majority of the recent high-lift research has been specifically aimed 
at understanding two-dimensional flows, both through experimentation and 
computation. The experimental studies have been particularly good at 
demonstrating the effects the geometry, Mach, and Reynolds numbers have on the 
airfoil lift and drag. Valarezo et al. has published a number of experimental 
investigations of three-element airfoils. His papers have been dedicated to slat and 
flap positioning [9,10] as well as some work on flow scale effects [11]. Valarezo's 
work has demonstrated adverse Reynolds number effects for multi-element airfoils, 
where the airfoil Clmax decreases as the Reynolds number increase. Chin et al. has 
provided some of the most detailed high-lift flow measurements in his 
measurements of the flow over a three-element airfoil [12]. This work has proved 
extremely useful to CFD code validation. An innovative, high-lift concept was 
investigated by Ross, Storms, and Carrannanto [13] and Storms and Ross [14] in 
their investigation of lift enhancing tabs on multi-element airfoils. Their preliminary 
findings suggest that in some cases Cl max can be improved with the addition of 
simple tabs near the trailing edge of an airfoil. 

There have been literally hundreds of attempts to analyze high-lift systems 
using coupled inviscid/ viscous methods-the few mentioned in this thesis are 
representative of most of this work. Most of the research used a two-dimensional 
potential method with an integral boundary layer scheme [15-17]. The primary 



advantage of this approach is the relatively low computational requirements. This 
makes these methods very attractive to designers who perform many iterations in 
the design process. More recently, quasi-three-dimensional versions of 
inviscid/viscous methods have been introduced [18]. This method combines a 
three-dimensional lifting surface method with two-dimensional viscous modeling to 
compute high-lift flows. The limit to this approach lies in the fact that the complex, 
three-dimensional viscous effects are represented with simple two-dimensional 
models. Also, with the three-dimensional potential code the computational 
requirements increase significantly. One of the more successful inviscid/viscous 
approaches is the MSES code [19]. MSES solves the two-dimensional Euler 
equations simultaneously with integral boundary layer equations. This method has 
given very good results, but is again limited to two dimensions. 

With the recent availability of quality experimental data, a great deal of CFD 
code validation work on multi-element airfoils has been performed. Due to the 
inherently viscous nature of the high-lift flowfield, most of the recent numerical 
investigations have solved some form of the Navier-Stokes equations. Perhaps the 
most widely used approach involves the solution of the incompressible Navier- 
Stokes equations on structured. Chimera grids [20-23]. Others have utilized the 
flexibility of the Chimera grid scheme to solve compressible, full [24] or thin-layer 
[25, 26] Navier-Stokes equations. Recently, some investigators have successfully 
solved Euler [27] and compressible Navier-Stokes equations on unstructured meshes 
[28]. All of the Navier-Stokes solutions mentioned here have demonstrated some 
ability to accurately resolve high-lift flows. The vast majority of the two- 
dimensional computational research discussed here was dedicated to analysis of 
conventional high-lift airfoils. A notable exception to this is Carrannanto et al/s 
analysis of the lift-enhancing tab concept [23] experimentally investigated by Ross, 
Storms, and Carrannanto and Storms and Ross [13, 14]. 



There are a few papers that do not fit neatly into one of the above categories, 
such as the papers that focus on specific flow features. For example. Savory et al. 
[29] discusses flow physics in the gaps between elements of a high-lift system,' 
Horton [30] dedicates his paper to the physics of high-lift flow separation, and van 
Dam et al. [31] discusses boundary layer transition and relaminarization on high-lift 
wings. Descriptions of the respective phenomena are included along with 
speculations about the driving mechanisms. As the understanding of high-lift flows 
increases, more of this type of research is sure to appear. Another paper topic is 
based on the development or assessment of high-lift system analysis techniques [32]. 
In a different type of paper, a summary of high-lift system design at Deutsche- 
Airbus is given by Flaig and Hilbig [4]. 

The lack of three-dimensional research becomes obvious when surveying the 
literature. This is made even more surprising by the fact that most of the survey 
papers discussed above indicate three-dimensional flow effects as a major driving 
force in high-lift system performance. Many of the survey papers contain some 
three-dimensional test results, but these offer only integrated effects without the 
details needed to begin to analyze the flow physics. On the computational side of 
things, very little has been done in the way of modeling three-dimensional high-lift 
systems. Most recently, Rogers [21] performed a preliminary Navier-Stokes analysis 
of a T-39 Sabreliner wing equipped with a high-lift system. Studies that 
computationally investigated complete aircraft with flaps deployed have lacked the 
grid resolution to adequately resolve the flow about the high-lift systems. They do, 
however, offer hope that a high-lift flowfield could be analyzed computationally 
with the proper grid resolution. 



1.3 SUMMARY OF THE PRESENT STUDY 


A CFD approach was chosen for this study because it may provide detailed 
flow information. Due to the viscous nature of high-lift flows, the simplest 
governing equations that could be used were the incompressible Navier-Stokes 
equations. It has been shown that compressibility can become an issue in high-lift 
flows for freestream Mach numbers as low as 0.15 [7], but for this study an 
incompressible solver was used due to the large savings of computer time. For high- 
lift flows at lower Mach numbers (Moo £ 0.2), any appearance of compressibility 
usually occurs in the vicinity of the slat element. The configurations investigated 
here were all two-element wing sections without slats. Therefore, it was felt that any 
loss of accuracy involved with assuming incompressible flow would be kept to a 
minimum. 

The goal of this study was to apply current computational tools to a simple, 
three-dimensional high-lift system. The work began with a two-dimensional 
computational study of the flapped airfoil that represented the basic geometry for 
the entire study. The two-dimensional work allowed for a grid resolution study and 
established a baseline set of results with which to compare at later stages in the 
investigation. The first three-dimensional computations were performed for a wing 
that fully spanned the test section of a wind tunnel. This was essentially the same 
configuration as the airfoils studied in two dimension, except the three-dimensional 
flow solver was used. These results were compared to the two-dimensional 
computations and experimental data to validate the approach. The major focus of 
the research was the computation of the flow over a wing with a half-span flap. A 
grid scheme was developed to adequately resolve this flow, so the details of the flow 
physics could be investigated. 

This project accompanies a wind tunnel test of the same geometry scheduled 
to occur in the NASA Ames Research Center 7-by-10 Foot Wind Tunnel. The CFD 



was done prior to the experiment so that the computed results could be considered 
when selecting a run schedule and locating the pressure sensors on the model. 
Every attempt is made to match the experimental set-up in the computations so that 
a meaningful comparison of the results can be made. This includes matching the 
flow conditions and geometries as closely as possible. For this reason, the wind 
tunnel walls were modeled in all of the computations. False walls will be needed in 
the experiment to shield instrumentation and mounting hardware from the flow, 
effectively reducing the size of the test section to 5x10 feet [14]. This 5x10 foot test 
section is the one modeled in the computations. The reference chord of the physical 
wind tunnel model (unflapped section) is 2.5 feet. The flow at the test section of the 
tunnel is at a freestream Mach number of Moo = 0.2 and a Reynolds number, based 
on the unflapped chord, of Rec = 3.7X10^. All length dimensions in the 
computations are non-dimensionalized by the model chord. 

This thesis contains a description of the geometry, governing equations, flow 
solver, turbulence model, and boundary conditions used for this study. The grid 
strategy developed to resolve flow about a three-dimensional high-lift system in an 
efficient manner is discussed. Computations of flow over a two-element wing 
between wind tunnel walls are compared with experimental results to verify the 
computational approach. The results for a half-span flap wing are also presented. 
Conclusions about the computational approach and flow physics are made. 



CHAPTER 2 


GOVERNING EQUATIONS 


2.1 INTRODUCTION 

This chapter covers the theoretical aspects of the research. The governing 
flow equations and related assumptions are discussed in Section 2.2. The complete 
derivations are left out of this chapter and can be found in Reference [33]. Section 
2.3 is dedicated to a non-dimensionalization and transformation of the governing 
equations which make the equations easier to solve numerically. A description of 
the Baldwin-Barth turbulence model makes up the final section of this chapter. 



This study involved solving the three-dimensional, incompressible Navier- 
Stokes (INS) equations. These equations make up a set of four partial differential 
equations which are derived from the laws of conservation of mass and momentum . 


In tensor form, the equations can be written as follows: 


— = 0 
dxi 

dui d(iliUj) - 6ij dp ,dw duj . 

+ — — + v( + — -) 

dt dXj p dXi dXj dxi 


(i) 


where i,j =1, 2, 3. The parameters in Equation (1) are defined as follows: 

xi represents the three Cartesian coordinates (x,y,z), 
w represents the velocity components (h,v,h>), 
p is the flow static pressure, 
t is time. 
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p is the flow density, 

v is the kinematic molecular viscosity of the fluid, and 
6ij is the Kronecker delta function. 

As shown. Equation (1) is classified as a mixed elliptic-parabolic set of equations 
[34]. 

The INS equations were selected because they are the simplest set of 
equations that could govern a high-lift flowfield. With the possible existence of flow 
features such as flow separation, viscous wake interaction, and confluent wakes and 
boundary layers, viscosity must be included in the modeling of the flow. Simpler 
forms of the Navier-Stokes equations (i.e. thin-layer or boundary layer equations) 
neglect aspects of the flow that were important to this study. It is not as obvious, 
however, whether the compressible or incompressible form of the Navier-Stokes 
equations should be used. It has been shown that compressible effects can appear in 
high-lift flows for freestream Mach numbers as low as Moo = 0.15 [7]. The 
compressibility appearing at such low Mach numbers is usually restricted to the slat 
element because the flow velocities are highest in this region. The simple geometry 
investigated in this study contains only main wing and flap elements, therefore 
reducing the possible effects of compressibility. Because the incompressible 
equations have been successfully applied to high-lift flows in two dimensions [20- 
23], it was decided that any loss in accuracy due to the choice of governing equations 
was minor when compared to the computational efficiency gained by using an 
incompressible flow solver. 

In general, there exists no closed form solution to the incompressible Navier- 
Stokes equations. Therefore, numerical methods are required to obtain a solution. A 
finite difference approach is used in this study. Finite difference methods replace 
the flowfield with a finite number of points (grid points) and solve the equations at 
these points. The primary concern with this method is that the equations will only 
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resolve flow phenomena that have characteristic scales greater than the distance 
between grid points. For laminar flows, this is not a problem because the smallest 
scales are on the order of the boundary layer thickness. For turbulent flows, the 
problem becomes much more difficult. The smallest scales in turbulent flows, the 
Kolmogrov microscales, are far too small to be resolved with grid points for all but 
the simplest flows [34]. Not only would points be needed to resolve these scales in 
the boundary layer, as with the laminar case, but this resolution would be required 
at every place in the flow containing turbulence. There is disagreement about the 
exact spacing required to resolve turbulent flow features, but one estimate is that 10 5 
grid points would be required to resolve 1 cm 3 of a typical flow field [34]. Clearly 
this is not currently feasible for most practical aeronautical applications. Therefore, 
additional assumptions are needed. 

One way to overcome this difficulty is to assume that physical flow properties 
consist of two components, one component associated with the mean flow property 
and the other associated with a turbulent fluctuation. For example, the velocity 
would be expressed as 


u « u +u' 


(2) 


where u is the total velocity, u is the mean velocity, and u' is the unsteady turbulent 
velocity fluctuation. All of the parameters in Equation (1) are replaced with this 
representation, and the entire equation is time averaged. The resultant set of 
equations are the well known Reynolds Averaged Navier-Stokes (RANS) equations 
[34]: 


dm 

dxi 


— 0 


( 3 ) 



dm 

dt 


+ 


d(uiui) - 6ij dp 1 

m + — \Tij 

dXj p dxi p 


pu'iu'j) 


where nj = p(— + —) and ^ . 

\ dxj dxi) p 

A note concerning this time averaging is in order. The time over which to 
average the equations is chosen to be large enough to capture small scale turbulent 
effects but small enough not to average out large scale unsteadiness in the flow. In 
practice, this is not a problem because the small turbulent time scales are much 
shorter than those associated with the flow in general. From examination of 
Equation (3) it is seen that Reynolds averaging process changes the equations in two 
ways. First, all the parameters in Equation (1) are replaced with an averaged value 
indicated by the overbar. The second difference is the addition of an extra term 
( “ pu'iu'j ) in Equation (3) know as the Reynolds stress. This term is "stress-like" 
because it adds "viscosity" to the flow due to turbulent fluctuations. The Boussinesq 
approximation uses this stress-like behavior to form the following approximation of 
the Reynolds stress based on the laminar stress relation: 


— — - dm duj 

-pu iU j *= Vt{ + ) 

dxj dxi 


( 4 ) 


where v, represents the turbulent, or eddy, viscosity. When Equation (4) is 
substituted into Equation (3), the following relationship is obtained: 


du 

dt 


i d(muj) -dij dp , ,dm duj s 

- + — - — - — + 0 + vt )( — + — -) 

dxj p dxi dxj dxi J 


( 5 ) 


In Equation (5) all of the effects of flow turbulence are modeled by a term, vt, that 
acts like additional molecular viscosity. But vt is an additional unknown, requiring 



12 


an additional equation, known as a turbulence model, to make the set determinate 
once again. This is known as the closure problem. The turbulence model is 
discussed in Section 2.4. All of the equations hereafter are assumed to be time, 
averaged unless specifically stated otherwise, and the overbars will be omitted for 
convenience. 


2,3 TRANSFORMED EQUATIONS 

The Reynolds Averaged Navier-Stokes equations as presented in Equation (5) 
are in a form that can be solved using numerical techniques. In this form, though, 
flows over general bodies would be very difficult to obtain because the current 
coordinate system would be inconvenient for a practical CFD problem. There is a 
non-dimensionalization and a transformation that will greatly simplify the solution 
of the RANS equations for complex configurations. First, the equations will be non- 
dimensionalized, and second they will be transformed into a generalized set of 
curvilinear coordinates. 

Many possibilities exist for non-dimensionalizing the RANS equations. The 
approach presented here is based on the implementation in the INS3D-UP code [35]. 
The following expressions are substituted into Equation (5) 


Ui _ Xi 
Ui = — , Xi = — , t 
Uref Xref 


tUref ~ P - Pref 

> P * 2 ' 

Xref P U ref 


Tij V _ , 

Tij = — v » « Re 

P U ref XrefUref 


( 6 ) 


which give 



dXi 

dui [ d(wuj) - 6ij dp 

dt dxj p dxi 


din 

+ (u + w,)(— + 

dXj 



( 7 ) 



As with the time average notation, the tildes will be dropped hereafter, but the 
variables remain non-dimensionalized unless otherwise stated. As shown, the non- 
dimensionalization does not change the form of the equations, but does allow for 
more convenient scaling between problems. 

A particularly useful form of the equations can be reached by combining 
Equation (13) and Equation (14) along with some rearrangement. 

BU BE BF BG n 

— + — + — + — «= 0 ( 8 ) 

Bt Bx By Bz ' * 


where 



w 



6 — €v f 

F - 

f-fi. 

G = 

g-gv, 


U 


V 


W 


u 2 + p 


UV 


UW i 

e = 

UV 

/ /- 

1 

/ g = 




+ 

1 

vw 


UW 


vw 


* 

+ 

^3 

Lmm - 



‘O’ 


■0 


■O' 


Txx 

, fv * 

Xyx 


Tzr 

e v = 

T. xy 

Xyy 

, g v ~ 

T zy 



Xxz 


Xyz 


X zz 


( 9 ) 


where e, f, g are convective flux vectors and e v , fv, gv are the viscous flux vectors. 
This form of the equations is known as conservation law form. Four equations have 
been combined into a single vector equation making its use easier from the 
standpoint of derivation and implementation. 
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Because the RANS equations. Equation (8), are to be solved using finite 
difference techniques, the computational mesh has a significant influence on the 
solution scheme. For example, if the flow problem is discretized with a mesh that is 
homogeneous (equal grid spacing in all directions) the finite difference schemes 
simplify dramatically. However, the number of physical problems that lend 
themselves to this type of representation are limited. Most applied flow problems 
require meshes which are quite curved as shown in Figure 3. The complexity of 
finite differences increases so significantly with an irregular mesh that most flow 
solvers transform the curved mesh into a homogeneous one. The flow equations are 
transformed and solved on the new mesh in the computational domain, and the 
solution is inverse transformed to yield results in the physical domain. This 
transformation certainly complicates the solution scheme, but the increase in 
complexity is less of a penalty than solving the equations on an irregular mesh. The 
real difficulty comes from the fact that every mesh must be conformed to suit its 
particular problem eliminating the possibility of an analytic expression for the 
general transformation. Therefore a general transformation is assumed and 
evaluated numerically during the solution process. 

For a three dimensional flow, a transformation of the form 

£ = £(x,y,z), (10) 

v = n(x,y,z), 

£ - U x >y,z) 

can be assumed. In this case, t], and £ represent orthogonal axes in the 
computational domain. Also, the mesh is assumed not to vary with respect to time. 

A two dimensional relation between the computational and physical planes is 
illustrated in Figure 4. Before the equations can be applied to the computational 



space they must be expressed as functions of the new space variables. To do this, 
consider the chain rule 


_d_^ _d__ 

dx 

d t d 
dy y d% 
d d 

Tz~ h Hi 


a „ a 

+ r}x-— + &—, 

an <?£ 

d „ d 
dt] d£ 

d . a 

dr] v <?£ 


( 11 ) 


where abbreviated partial derivative notation is used (|x a — , etc.). Equation (11) 

dx 

can be substituted into the RANS equations to transform them into computational 
space. Before the real usefulness of the transformation can be appreciated, a little 
more needs to be said concerning the metrics of the transformation (&,rj y , etc.). 
Expressions for the metrics are needed in order to evaluate transformation. To 
derive expressions for the metrics, consider the following differential expressions 
[ 34 ]: 


dx = + x v dr] + xsd £ , (12) 

dy - ysd% + yr4t] + y J, 
dz » zi^|f + z n dr] + zzdt, 

and similarly 

d% = ^xdx + %ydy + l=zdz, ( 13 ) 

dr] ■= r]xdx + r]ydy + rjzdz, 
d£ = £xdx + £ydy + &dz 


which can be put in the matrix form 
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A = BC, (14) 

C *= DA 


if A, B, C, D are defined as follows: 



' dx 


Xr, X£ 

Am 

dy 

/ B- 

yz y* yt , 


dz 


zt Zr) zz 


d£ 


& 1 y & 

C = 

dr] 

/ Z) = 

r]x rjy 7] z 


dt, 


& & & 


(15) 


It follows that 


B' l A - * C 


(16) 


and since 


C~DA 


(14) 


the result 


B' 1 = £> 


( 17 ) 


is attained. If B is inverted , expressions for the metrics are attained [34]. 
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D - B~ l «■ J 


y^zt - ytZn 
-(ysZi-ytZt) 
ygn - y& 


-(XrZt - XtZr,) 
X(Zt - x& 
-(XfZr, - X*Zt) 


x n y t - x t y n 
-(x t yt-x c y t ) 
x { y„ - x v y t 


J is the Jacobian of the transformation given by 


j = *?>£) = J d ( x >y> z ) 

d(x,y,z) / ^(f,i7,S) 


X$ 

x n 

XK 

ys 

y* 

yz 

Z| 

Zrj 

Z£ 


(18) 


(19) 


Obviously, the determinant in Equation (19) must be non-zero for the 
transformation to exist. Physically, this prevents the mapping of a zero or negative 
cell volume in one domain to a positive cell volume in the other. The individual 
metrics in Equation (15) can be solved if the expressions in Equation (18), or more 
specifically the partial derivatives x^y^zz, etc., are known. There are two ways of 
finding values for the partial derivatives. They can be computed analytically if a 
closed form transformation (or inverse transformation in this case) exists. In general, 
however, this relationship does not exist. The alternative method requires the use of 
finite differences to evaluate the derivatives. The finite difference approach is used 
here and is discussed in the section on differencing (Section 2.4.2). For the rest of 
this section it is assumed that all of the metrics are known. With known metrics. 
Equation (8) can be transformed by using Equations (29-31). 


BU BE BE . dE _ BF 

+ Tlx + Lx 1- |fy 

dt 3% dt] v d% * dt 

3F BF dG dG . dG „ 

+ T]y h 1- Tlz h Lz ■« 0 

Brj B£ B% Bt) B£ 


( 20 ) 
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This equation can be put back into conservation law form by employing th,e 
following definitions [36]: 


Ei = j(E& + F§ y + Gb), 
1 

Fi m —(Et]x + Frjy + Grjz), 
Jf 

Gi = j(E& + F£, + G&) 


( 21 ) 


Equation (20) becomes 


dU i dE i dF\ 

+ + + 

dt dr) 


dGi 


~ 0 


( 22 ) 


Recall that the vectors E, F, and G were made up of a convective flux and a 
viscous flux. The viscous terms contain partial derivatives that must also be 
transformed using Equation (11). The transformation of the viscous fluxes is 
omitted here but the details can be found in Rogers [35]. The following is a 
summary of the results: 


I*-! 

J 


O' 

u 

v 

w 


J J J 


%xU + %yV + % Z W 
Z-x(p + U Z ) + gyllV + |zWW 
gxUV + %y(p + V 2 ) + %zVW 
%xUW + + %z(p + W 2 ) 


evi 


0 

^xTxx + ifyTxy + l^zTxz 
^xTxy + ^yTyy + ^zTyz 
%xTxz + yTyz + ^zTzz 


( 23 ) 
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/i 


Tjxll + r]yV + 7]zW 
T]x(p + U 2 )+ T]yUV + YjzUW 
r\xUV + rj y(p + V 2 ) + TJzVW 
Tjxuw + rjyvw + rjz(p + w 2 ) 
£rW + V + £zW 
£x(p + U 2 ) + £yWV + t^zUW 
&MV + g>(/7 + V 2 ) + &VW 
&cUW + £yVW + £z(p + W 2 )\ 

0 du du du . 

+r7x— + ^— ), 
o>?7 <?£ 

_ 0 dv dv , <?V 

r>y - 2 u(^— + T]y — + gy — ), 

dt, drj dt, 

~ n rf- dw dw dw. 

rzz = 2 u(|z— + r/z— + &— ), 

d% dr] d£ 


gi = 


Txx 


f* 


gvl 


0 

rjxTxx + VjyTxy + TTJzXxz 
rjxTxy + TTfyXyy + TJzTyz 
TJxTxz + rjyTyz + TjzTzz 
0 

f^xtxx + t^yXxy + t^zXxz 
£xTxy + £yXyy + £zXyz 
£ xXxz + £yXyz + £zXzz 


dw 


r/Jr ^ du - du s . - i /rr 

- + T,x ^ + ^ ) + + 

du du du dv 

Tw = vK ^ + + & + ®ig + 

* " [(?y ^ + *7^ + ^ +(M, li + 


dw dw s 

T]z +& )], 

dr] ^dS 
dv dv 

n, ^ + ^ )h 
dw dw 

dr] * d£ n 


2.4 B ALDWIN-BARTH TURBULENCE MODEL 

The one-equation Baldwin-Barth turbulence model [37] was used for all of the 
cases discussed in this thesis. A number of researchers [22-24] have shown that this 
model does a reasonable job in predicting high-lift, multi-element airfoil flows. 
Specifically, Mani and Bush [24] have demonstrated a dramatic improvement in 
results computed with this model versus those performed using the algebraic 
Baldwin-Lomax model [38]. Rogers et al. recently studied four one- and two- 
equation models for use in computing multi-element airfoil flows [22]. He 
concluded that the all of the models tested "produced very similar results in most 
cases." This model requires approximately the same grid resolution as the algebraic 
models, which is significantly less than the multi-equation models require [37]. 
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These facts, along with the relatively low computational requirements, made this 
model very attractive for the use in this study. 

The current implementation of turbulence models in INS3D-UP requires that 
the code be run in a fully laminar or fully turbulent mode. All of the current cases 
were run in a fully turbulent mode, with no modeling of transition. This did not 
hurt the comparison with the experimental data too badly because transition strips 
were placed on the upper (x/c = .10) and lower (x/c = 0.05) surfaces of the main 
element in the experiments. No transition strips were used on the flap, but the large 
adverse pressure gradients typically seen on the upper surface of the flap near the 
leading edge are likely to cause natural transition close to the leading edge. Laminar 
flow is more apt to exist on the lower surface of the flap which would create an 
inconsistency between the computations and the experiment. This difference should 
have a minor effect on the flowfield as a whole. Transition strips will be used in the 
planned half-span flap experiment. 

The Baldwin-Barth model was developed through a simplification of the 
standard k~£ model equations [39]. The k-e model is a system of two partial 
differential equations for k, the turbulent kinetic energy production, and e, the 
energy dissipation rate. To combine the k - e equations, a new parameter, Rt, is 
defined that includes the turbulent kinetic energy production and dissipation terms, 

k 2 

Rt “ — (24) 

V£ 

where Rt is known as the turbulent Reynolds number. With this definition, the two 
partial differential equations were reduced to a single partial differential equation 
for Rt. The equation is iteratively solved for the entire flowfield using a line 
relaxation scheme. Once Rt is known, the eddy viscosity, vt , can be computed from 
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v t - C(vRT)Dfi 2 (25) , 

A-l-exp(-//A + ) 

A-l-expC-//^) 

y + =u x y I v 

where w T is the friction velocity ^jr wall / p walI/ C = 0.09, = 26, and A> + = 10. C, 4 + , 

and are empirically determined constants [37]. 



CHAPTER 3 


NUMERICAL METHODS 

3.1 INTRODUCTION 

All of the computed results in this study were attained using the INS2D-UP 
and INS3D-UP (Incompressible Navier-Stokes with Upwind Differencing) flow 
solvers [35]. The INS codes provide a means for solving the Incompressible Navier- 
Stokes equations in a manner that is fast, robust, and easy to perform. This section is 
dedicated to describing the numerical solution scheme incorporated in INS2D-UP 
and INS3D-UP. The three-dimensional method is discussed because the two- 
dimensional method is a subset of this case. Only the aspects of the code pertinent 
to this research are discussed. 

3.2 ARTIFICIAL COMPRESSIBILITY FOR STEADY-STATE PROBLEMS 

INS3D-UP solves the Incompressible RANS equations in their primitive 
variable form. To enhance convergence, the method of artificial compressibility is 
employed [40]. This method simulates a compressible flowfield by introducing 
artificial pressure waves which propagate through space. The implementation of 
this method involves adding an artificial compressibility term to the continuity 
equation. 


^ = -{SVV (26) 

In this case, (3 represents the speed at which the pressure waves travel. It is stressed 
that this is not a physical phenomena, but merely a numerical scheme used to 
enhance convergence. This relationship does not affect the final solution because the 


.22 
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pressure waves (~) disappear as the solution approaches steady state. As a result, 

the steady continuity equation remains unchanged. For a steady solution, the time 
variable t does not represent physical time, so it can be replaced with a pseudo time 
term r [35]. This is done to avoid confusion in cases where a time-accurate solution 
is needed. With the addition of (26) to the terms in the RANS equations (22) become 


dU i dEi dFi dG\ 

+ + + 


drj dt, 


ei — evi 


Z 1 - Q gl - gvl 

J ' ^ J ‘ 


P(^xU + §yV + %zW) 
gx(p + U 2 )+ gyUV + %zllW 
%xUV + %y(j? + V 2 ) + gzVW 
gxltW + gyVW + %z(j) + W 2 ) 
P(T]xU + TjyV + t]zW) 

t]x(p + u 2 ) + rjyuv + rjziiw 
rjxuv + r] y (p + v 2 ) + rjzvw 

TJxUW + T]yVW + 7Jz(p + w 2 ) 
P(£xU + t,yV + £zw) 

£x(p + U 2 ) + £yUV + £zUW 
&UV + £y(/? + V 2 ) + frVW 
tpUW + £yVW + £z(p + W 2 ) 


^xTxx + ^yXxy + ^zXxz 
^xTxy + %yXyy + ^zXyz 
^xXxz + ^yXyz + ^zXzz 

r o 

TJxXxx + TjyXxy + TjzXxz 
TJxXxy + TjyXyy + V)zXyz 
TjxXxz + TjyXyz + TjzXzz 

^ 0 
t^cXxx + £yXxy + £zXxz 
^xXxy + £yXyy + ^zXyz 
^xXxz + £ yXyz + ^zXzz 


There are physical and numerical advantages to using the method of artificial 
compressibility [35]. Physically, the pressure becomes directly coupled to the 
governing equations through the artificial compressibility term (26). Most primitive 
variable solution schemes require the solution of Poisson's equation for pressure at 
the end of each iteration, resulting in an indirect coupling of pressure to the rest of 
the flowfield [34]. The direct relationship established by the artificial compressibility 
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term is not only more physically realistic, but it also speeds up convergence because 
Poisson's equation no longer needs to be solved each iteration. This dramatically 
reduces CPU time per iteration. In addition to speeding up convergence, the 
addition of artificial compressibility also changes the form of the equations from a 
mixed elliptic-parabolic set to one of a pure hyperbolic nature [35]. The hyperbolic 
classification opens up the possibility of upwind differencing and the use of 
marching schemes to solve the equations. These details are the subject of the 
following sections. 

3.3 DIFFERENCING 

Before the flow equations can be solved numerically, the partial derivatives 
must be approximated in some manner. One widely used representation of the 
partials in the aerospace field is the finite difference approximation. Such an 
approximation does two things: the flowfield is transformed from a continuum to a 
discrete set of points, and the partial derivatives are replaced by algebraic 
expressions based on Taylor Series expansions. These two results allow 
approximate flow solutions to be obtained that would otherwise be impossible. The 
differencing schemes used in the INS3D-UP code are discussed here. 

3.3.1 Transformation Metrics 

As previously discussed, analytic expressions for the transformation metrics 
(10) are generally not known. Therefore some approximations must be made. As 
with the governing equations, it is convenient to compute the metrics using finite 
difference expressions. This allows the choice of any right-handed, body-fitted 
coordinate system in the physical domain as well as a uniform, orthogonal grid in 
the computational domain. For a grid that does not vary with time, a second-order 



central difference can be used to compute the partial derivatives in the metrics. 
Equation (28) is such an example for a uniform, orthogonal grid. 


dx i,j,k 


as 2A£ 


( X i+l,j,k X i-\,j,k) 


(28) 


The individual partials can then be spatially averaged to form the metrics. For 
example, if an average operator /z is defined as 


x ij,k 


2AS 


( X i+l,j,k X i-l,j,k) 


(29) 


then the final form of the metrics becomes [35] 


& - ^(y?W« { ) - ^(y^(zrj) 


(30) 


3.3.2 Pseudo-Time Derivatives 

The flow equations are solved by marching in pseudo-time until a steady- 
state solution is reached. Therefore a representation of the pseudo-time derivatives 
is required. The flow equations have been put in the following form (Section 3.2): 


am 

dx 


BE\ dFi 

+ + + 

dS Brj 


BGi 


*0 


(31) 


or 


au i 

dx 


S3 —R 


(32) 


where 
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R = 


dE i dFi dGi 
d% dr] d£ 


(33) 


The pseudo-time derivative can be replaced with a simple, implicit Euler difference 
since accuracy in pseudo-time is not required [35]. 


a u ; +1 

At 


~-R 


n + 1 


A U? +l = U? +1 -U? 


(34) 


where n represents the pseudo-time level. The right hand side of (34) is then 
linearized to give 


I 

JAr 


(dR\ 

\dU) 




-R n 


(35) 


where I is a 4x4 identity matrix. The details regarding the formation of R n and the 
solution of the equations are left to subsequent sections. 


3.3.3 Convective Flux Differencing 

In the previous section, the residual term, R n , was defined as 





(36) 


where 


Ei = 


ei - evl Fi m /* - fa Gi = gl - gv 1 


J 


J 


J 


(23) 
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and 


SE X 1 de 1 


if V, 


■),etc. 


(37) 


This section will only cover the convective terms while the viscous terms are 
discussed in the following section. The convective terms are represented using an 
upwind difference formulation. This scheme expedites convergence and also 
increases the robustness of the code [35]. CFD codes that utilize central differencing 
typically require the addition of numerical dissipation to damp out numerical 
oscillations caused by the non-linear convective fluxes. This has been a problem 
facing the CFD community for years since the amount of dissipation added can 
influence the final solution. Upwind differencing, however, provides implicit 
dissipation to the equations, removing the need for any explicit, user-supplied 
dissipation. Also, an upwind scheme contributes to the diagonal term of the 
Jacobian of the residual, making the scheme nearly diagonally dominant [35]. A 
diagonally dominant matrix is much easier to numerically solve than a non- 
diagonally dominant matrix. The upwind scheme used in INS3D-UP employs flux 
difference splitting based on the sign of the eigenvalues of the convective flux 
Jacobian [35]. The details that follow are based on Roe's method for solving the 
compressible gas dynamic equations [41] and are taken from Rogers [35]. 

Each dimension of the problem is considered separately as follows: 


^ g l m e i+ 1/2 g /-l/2 


(38) 


where e is a numerical flux given by 





(39) 


E i( U i) is the value of E 1 obtained using the values of the primitive variables at U r 
0«+i/2 is a dissipation term. A first order upwind scheme results if <j>. +1/2 is defined 
such that 


0/+ 1/2 = i ^7+l/2 — ^'i+1/2 (40) 

where AZs* is the flux difference across ± traveling waves. The flux difference can 
be evaluated if the following definitions are established 

U " + U t ), U Mn = U M - U, (41) 

and if the similarity transformation is considered. The Jacobian can be expressed as 

A=XAX~ l (42) 

if A is the matrix with the eigenvalues of A on the main diagonal and zeros off the 
diagonal, and X is the matrix of eigenvectors of A. Since the eigenvalues of A 
correspond to the slopes of the characteristics, the direction of information 
propagation is determined by the signs of the eigenvalues. To model this 
phenomena in a physically meaningful fashion, it is necessary to distinguish 
between positively and negatively sloped characteristics (this is also necessary for 
numerical stability). Therefore, consider the following splitting of the eigenvalue 


matrix: 
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A + - XA + X ~ l , A~ - XA~X' 1 (43) 

where A ± is the matrix containing the ± eigenvalues. This splitting allows 
appropriate differencing based on the local direction of information propagation. 
Now the flux difference can be computed by 

M? +V2 -A*(U)AU i+l ,2 (44) 

Examples of higher order schemes, as well as the complete eigenstructure of A, can 
be found in Rogers [35]. 

3.3.4 Viscous Flux Differencing 

The final components of the residual are the viscous flux terms. These terms 
are differenced using second order central differencing in the INS3D-UP code. For 
example 


e vu+l/2,j,k e vii-l/2,j,k 




Af 


(45) 


where the values at the half-points are obtained by averaging neighboring values. 


e i + V2 “ “0/ + l + *f) 


(46) 


3.4 IMPLICIT SCHEME 

The flow Equation (35) is numerically solved using a block Gauss-Seidel 
iterative scheme with line-relaxation. The solver sweeps through the computational 
domain with varying direction. All of the terms along the current line are solved for 
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implicitly, while the terms not on the sweep line are multiplied by the current AC/" 
and moved to the right hand side. The direction of the sweeps alternates as 
prescribed by the user until the residual and the divergence of the flowfield 
approach zero. This solution strategy follows that outlined by MacCormack [42]. 

All of the information in Equation (35) is known except for the A C/" +1 term. 
The Jacobian of the residual vector is computationally very expensive to form, so it is 
approximated as described by Rogers [35]. At each iteration the AC/" +1 term is 
solved implicitly. When AC/" +1 is known, the flow variables can be updated using 
Equation (34). 



CHAPTER 4 

COMPUTATIONAL GRIDS AND BOUNDARY CONDITIONS 

4.1 INTRODUCTION 

This chapter covers the information relevant to the building of the 
computational meshes. Before the grids are discussed, a description of the 
geometries that were studied is included. An overview of the general grid building 
scheme, as well as a description of the grid building tools, makes up the third section 
of this chapter. Section 4.3 contains the information common to all of the grids that 
were built. Specific details for each grid are discussed in the sections dedicated to 
airfoil, full-span, and half-span grids. The final section identifies the different 
boundary conditions used in this simulation. 

4.2 GEOMETRY 

The geometries studied were all based on the NACA 632-215 Mod. B airfoil 
section [43]. This airfoil was fitted with two different riggings of a 30% Fowler flap. 
Figures 5-7 show the airfoil section with and without the flap elements. The two flap 
riggings studied were representative of a take-off and landing configuration. The 
rigging details are contained in Table 1. Both configurations were set at 10 degrees 
angle of attack. 


Condition 

Deflection (deg.) 

Gap /Chord 

Overlap /Chord 

Take-off 


0.030 

0.042 

Landing 

30.0 

0.023 

-0.0039 


Table 1, Flap Rigging Parameters. 

The gap and overlap are defined in Figure 8 and are expressed in units of unflapped 
airfoil chord. As shown, the gap is defined as the vertical distance between the 
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lower surface of the trailing edge and the highest point on the deflected flap 
element. Overlap is the horizontal distance between the main element trailing edge 
and the leading edge of the deflected flap. It is defined so that in a configuration 
with negative overlap the flap is positioned completely downstream of the main 
element. 

The full-span flap cases are a simple extrusion of the airfoil into a wing for 
which all of the spanwise stations are identical. The wing fully spans a modeled 
wind tunnel test section. The test section is modeled after the NASA Ames 7 -by- 10 
Foot Wind Tunnel with the false wall discussed by Storms and Ross [14]. The half- 
span flap cases consist of the flapped airfoil over half of the tunnel span with the 
unflapped section composing the second half of the span. The half-span flap 
geometry is shown in Figure 9. 

4.3 GRIP BUILDING PROCESS 

The grid building process begins by creating a surface grid from a surface 
definition. The difference between the surface definition and grid is largely a matter 
of detail. A surface definition consists of enough information (x,y,z locations) to 
completely and accurately define the surface shape, whereas a surface grid contains 
specific refinement needed to resolve the desired flow features. Also, in the case of a 
C-gnd, as used for the wing section grids, a wake cut is part of the surface grid. All 
of the surface grids in this study were created in two-dimensions and extruded into 
three-dimensions where necessary. 

The next step is to create a volume grid. The hyperbolic grid generator 
HYPGEN [44] was used for the grids in this study. HYPGEN begins at the surface 
and creates the volume grid by marching outward a prescribed distance using the 
solution of hyperbolic partial differential equations. The solution to these partial 
differential equations must satisfy two orthogonality relationships and one cell 



volume constraint [44]. A volume grid is generated in this fashion for each airfoil 
element. An example of this can be seen in Figure 10. 

The final overset grid was created using the PEGSUS code [45]. PEGSUS is 
based on the Chimera scheme of Benek, Buning, and Steger [46]. The Chimera 
scheme provides one method of joining individually generated grid zones into a 
single region. As part of the merging process, communication must be established 
between the zones so that each grid can be influenced by neighboring regions. Also, 
grid points in overlapping regions that do not make physical sense need to be 
removed from the computational domain. For example, grid points in one zone 
should not fall inside a solid surface of another zone. This is illustrated in Figure 11. 

In the Chimera approach, the grid zones communicate with one another 
through an explicit interpolation scheme. This interpolation occurs at grid 
boundaries, imposing a type of boundary condition. Each boundary point is 
updated by interpolating flow information from a surrounding cell in an 
overlapping zone. In Figure 12, point A is updated by interpolation flow 
information from the points marked B. At least one cell of overlap is needed for this 
interpolation link to be established by PEGSUS. The flow variables at the points 
requiring interpolation are updated after each flow solver iteration. This approach 
ensures a smooth solution across the boundaries, but does slow convergence due to 
changing values at the boundaries. 

PEGSUS allows the user to "cut holes" in a mesh to prevent non-physical 
situations as illustrated in Figure 11. When a hole is "cut" all of the points falling 
within the hole boundaries are removed from the computational domain. Each grid 
point contains an integer flag, in addition to the coordinates of the point, known as 
the Iblank value. For a normal point inside the computational domain, Iblank is set 
to one. When a point is removed as part of a hole, its Iblank value is set to zero (the 
point is said to have been blanked). A zero Iblank value flags the flow solver to 
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ignore this point when looping through the grid. An example of how holes are used 
to avoid problems in overlap regions is shown in Figure 13. The points that were 
previously inside the solid boundary have now been blanked and no longer appear 
in the computational domain. Since the grid in which the hole was cut is no longer 
continuous. Figure 14, additional information must be provided at the boundary of 
the hole. These points are known as fringe points. The flow variables at the fringe 
points are found by interpolating from other grids in the manner described above. 
Care must be taken when defining the holes so that at least one cell of overlap exists 
at all hole boundaries. In Figure 12 sufficient overlap exists, and point A will be 
updated by interpolating the values at the four B points. 

The grids are merged and interpolation stencils are formed by the PEGSUS 
code. The user provides the individual grids and a PEGSUS input file. The input 
file contains the definition of hole and outer boundaries in addition to zone location 
information. Using the zone positioning capabilities of PEGSUS, individual zones 
can be translated and rotated so that they are correctly placed with respect to the 
other zones. In this study, the model angle of attack was set using this feature. An 
interpolation database is created for all of the points that require updating by 
interpolation and is written to a file. The multi-grid and interpolation database files 
are used by the flow solver. 

4.4 AIRFOIL GRIDS 

The grids used in this study were modeled after those used by Carrannanto et 
al. in his study of multi-element airfoils [23]. In this study, however, wind tunnel 
walls were also modeled in the computations. This was done to obtain results that 
would be more directly comparable with wind test results. C-grids were used to 
model the airfoil elements and an H-grid was used to model the wind tunnel walls. 
A grid resolution study was performed to investigate the effects of grid density on 



the flow solution and to establish grid density requirements for the three- 
dimensional work. It was found that the resolution used by Carrannanto et al. was 
adequate to resolve the flow physics up to angles of attack near Clmax- Since it was 
desired to have the same cross-sectional grid for the two- and three-dimensional 
cases, and since the grid size must be kept as small as practical for three-dimensional 
studies, no grid points were added. 

The grid sizes for the main and flap elements were 227x80 and 155x59 
respectively. The circumferential grid spacing was 10 -3 chords at the trailing edge of 
both elements. The wall spacing was 10“ 5 chords (y + ave~l) for both of these zones. 
In the normal direction, the main element grid extended approximately two chords 
while the flap grid was on the order of 0.15 chords. Wake cuts extended 2.5 chords 
downstream for the main element and 0.3 chords for the flap. The wake cuts for 
both elements were located in the manner suggested by Carrannanto et al. The 
wake of the main element followed the upper surface of the flap and was offset by 
the gap between the elements (see Figure 15). Once aft of the flap element, the wake 
cut made a straight line parallel to the chord of the deflected flap. The wake cut of 
the flap grid extended from the trailing edge of the flap along the flap chord line. 
Elliptic smoothing was applied to the flap grid aft of the flap trailing edge as 
suggested by Carrannanto et al. Smoothing expands the circumferential grid lines in 
the vicinity of the wake cut, reducing the aspect ratio of the grid cells in this region. 
The result is accelerated convergence and improved interpolation between grid 
zones [23]. Trailing edge thickness was neglected in the computations for both 
elements. It was felt that any effects of the blunt trailing edges would be local in 
nature and would not reduce the overall accuracy of the computations. More recent 
results have confirmed this assumption [47]. 

The wind tunnel grid was designed for computing the effects of inviscid wind 
tunnel walls on the model. The size of the grid was 85X30, and resolution studies 



showed very little sensitivity to the dimensions of this grid. The grid extended ±15 
chords in the streamwise direction. The walls were located two chords above and 
below the center of the airfoil at zero degrees angle of attack. The vertical distance 
between the walls was constant along the length of the tunnel grid. The physical 
wind tunnel walls diverge slightly through the test section to account for the growth 
of the boundary layer, and this divergence was neglected in the computations. 

The three zones were combined into a two-dimensional multi-zone using 
PEGSUS. The complete grid contained approximately 30,000 grid points. An 
example of an airfoil grid can be seen in Figure 16. This grid represents the landing 
configuration at ten degrees angle of attack. 

4.5 FULL SPAN FLAP GRIPS 

The full-span flap grids provided a direct extension of the airfoil work info 
three-dimensions. The two-dimensional planes were copied and stacked together in 
a spanwise row to form a three-dimensional grid prior to the merging of the zones 
using PEGSUS. This allowed for different spanwise resolution in each zone. The 
work in two-dimensions allowed all of the Chimera details to be worked out prior to 
creating the three-dimensional grid. This offered significant savings of both time 
and computer resources. Part of the reason for studying this configuration was to 
determine the number of spanwise planes needed to resolve any three dimensional 
flow. As the study progressed, no strong effect of spanwise resolution was seen for 
this case. The final grid sizes were 227x30x80, 155x30x59, and 85x30x30, for a total of 
900,000 points. The physical spanwise dimension was determined by the wind 
tunnel size, including the presence of false walls needed for the experiment [14]. 
The span of the grid was 2 chords. 
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4.6 HALF SPAN FLAP GRIDS 

The half-span grid was composed of six zones. Cross sections of three of 
these zones, the flapped main element, flap element, and wind tunnel grid, were 
identical to those used in the two-dimensional computations. The flapped part of 
the wing extended from the mid-span of the tunnel to a wind tunnel wall. The other 
half of the tunnel was spanned by the unflapped wing element. This cross section 
was modeled using a C-grid of dimensions 185X52. As with the other wing sections, 
the wall spacing was 10*5 chords with a circumferential trailing edge spacing of 10"3 
chords. The wake cut of this grid extended 2.5 chords behind the wing in a direction 
parallel to the section chord. The outer boundaries were located a normal distance 
of two chords from the surface. 

The remaining two zones were needed to model additional solid boundaries. 
One such surface was the exposed region where the flapped main element met the 
unflapped main element at the mid-span. The flapped section has a smaller chord 
than the unflapped section as well as a cut-out, or cove, region where the flap is 
stowed when the high-lift system is not in use. When the flap is deployed, this 
mismatch area is exposed and must be modeled. Figure 17 illustrates the location of 
this surface. The grid used to represent this surface was an H-grid deformed to 
match the airfoil contours. This grid is hereafter referred to as the patch grid. As 
with all of the grids, the patch grid was initially generated as a two-dimensional 
plane. The plane contained 33x33 grid points. Figure 18 shows a patch grid plane. 
The final zone was needed to model the flap edge surface. The grid used to do this 
was also an H-gnd deformed to match the flap cross section. This plane contained 
36x40 points. A flap edge grid plane is shown in Figure 19. Both of these grids were 
expanded outward to establish an overlap region so that interpolation could occur. 

As with the full-span flap case, all of the grid planes were copied into 
spanwise stacks to fill out the third dimension. Forty spanwise grid planes were 
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used for the following grid zones: flapped and unflapped main elements, flap 
element, and wind tunnel. The patch and flap edge grid zones each contained ten 
spanwise planes. The three dimensional grids were then passed to PEGSUS for 
merging and establishing the interpolation links between the zones. Proper Chimera 
interpolation between grids proved a challenge in two locations; the overlap region 
between the flapped and unflapped main elements and the flap edge region. 

The mid-span location was complicated because two grids of different shapes 
were required to fit together and communicate. As mentioned, the shape difference 
between the main element sections was corrected through the use of the patch grid. 
The first spanwise plane of the patch grid was located at the solid surface, coincident 
with an unflapped grid plane, and the grid was marched 1.5X10"^ chords into the 
flapped portion of the grid. The interpolation between the main element grids was 
more difficult to establish. The unflapped main element extended 10 “ 3 chords and 
four grid planes into the flapped element grid. This overlap allowed the end planes 
of each grid to receive interpolated information from the interior of the other grid. 
This overlap created an additional problem, however. The inner plane of a C-grid is 
a computational boundary, hence requires boundary information for a solution to 
exist. For the most part, the flapped and unflapped surfaces were coincident, so 
solid surface boundary conditions could be used. Near the trailing edge, this 
approach did not work because a portion of the flowfield in one grid (flapped main 
element) would have coincided with a solid surface in the other grid (unflapped 
main element). This would have created an inconsistency in the solution. To 
circumvent this, the portion of the unflapped trailing edge that extended into the 
computational domain of the flapped grid was defined as a Chimera boundary, and 
thus was updated by interpolation. Figure 20 shows a cross section of the grid in the 
overlap region with the appropriate boundary conditions shown. 



The flap edge region was simpler to set up, but turned out to play an 
important role in the convergence of the solution. The first approach was to place 
the last grid plane of the flap edge grid at the solid surface and march the grid away 
from the surface as was done with the patch grid. Using this grid scheme, the 
solution would not converge. This was traced to an unsteady flow separation 
occurring at the comer of the flap edge and flap surface. The flap edge grid was not 
able to resolve this feature with only boundary information. The problem was fixed 
by extending the flap edge grid three grid planes into the flap surface. The points 
inside the surface were blanked from the computation and the solid surface 
boundary condition was applied to the plane that was third from the end of the grid. 
This allowed the flap edge grid to resolve the flow separation at the comer, and the 
solution converged. 

No grid resolution study was performed for this grid due to the 
computational requirements for a single solution. Because the solution captured 
small scale flow features, such as secondary separations at the flap edge, the 
resolution appeared sufficient for this study. The final grid contained 1.8 million 
grid points. Figure 21 shows the final surface grid for the half-span flap geometry in 
a landing configuration. 

4.7 BOUNDARY CONDITIONS 

All of the solid surface points on the wing were modeled using a no-slip 
boundary condition. This condition specifies zero velocity and a zero normal 
pressure gradient at the surface. Because a C-gnd topology was used, wake-cut 
boundary conditions were also needed. Wake points are updated by a first order 
averaging of the points on either side of the wake cut. The wind tunnel walls were 
represented by slip-walls which impose a zero normal gradient for all flow 
variables. The tunnel inflow condition was prescribed with uniform normal velocity 
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and constant total pressure. The outflow of the tunnel consists of a constant static 
pressure and extrapolated velocity. The outer boundaries of all grid zones falling 
inside the wind tunnel grid were Chimera boundaries. 



CHAPTER 5 


RESULTS AND DISCUSSION 

5.1 INTRODUCTION 

This chapter includes a presentation of the results of this research. The Airfoil 
Results section establishes a starting point for this study. Much initial work was 
done in two-dimensions due to the faster turn around time associated with this 
simplified problem. The full-span flap investigations provided a good validation 
case for the three-dimensional flow solver. This case also allowed for the evaluation 
of the grid development work performed in two-dimensions. The final section of 
this chapter provides a discussion of the interesting results seen in the half-span flap 
computations. Specific attention is given to the differences between the full- and 
half-span cases. 

5.2 TWO-DIMENSIONAL AIRFOIL RESULTS 

To make comparisons between solutions meaningful, it was important to 
establish a convergence criteria so that all solutions could be converged to the same 
level. The INS2D-UP codes writes the maximum divergence of velocity value at 
each iteration. From a physical point of view, the divergence represents how well 
mass is conserved for the incompressible flowfield. With the artificial 
compressibility term added to the continuity equation, the divergence values 
indicate the levels of unsteadiness in the solution. As the divergence approaches 
zero, the solution approaches steady-state. For this work, the solution was iterated 
until the divergence of velocity fell to the order of 0.1. This typically corresponded 
to a residual reduction of five orders of magnitude. Convergence was reached in 
300-500 iterations for angles of attack below Q max - Near stall, up to 1200 iterations 
were sometimes required to obtain a steady-state solution due to the inherent 
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unsteadiness in the flow. No computations were performed beyond stall. A 500 
iteration run using 30,000 grid points required approximately eight minutes of CPU 
time on a Cray C-90 supercomputer. 

The computed two-dimensional results are compared to data taken in an 
unpublished study at NASA Ames Research Center. The experiment occurred in the 
U. S. Army 7- by 10- Foot Wind Tunnel. The flow conditions for the test were 
freestream Moo = 0.2 and Re = 3.7xl0 6 based on the undeflected airfoil chord. The 
model had an unflapped chord of 2.5 feet. Due to mounting hardware, false walls 
were placed in the tunnel, reducing the effective test section span from seven to five 
feet. The lift coefficients were computed from an integration of pressures measured 
at the model mid-span. Boundary layer control was used at both walls to maintain 
two-dimensional flow at the wing-wall junction. Three blowing slots were located 
on each side wall, two above the upper surface of the main element and one above 
the upper surface of the flap. Air was injected along the wall parallel to the wing 
surface. The blowing rate was adjusted until the pressures at three spanwise 
stations agreed. The data used for comparison was not corrected for any wind 
tunnel wall effects. 

Figures 22 shows the variation of lift coefficient with angle of attack for the 
10 flap deflection (take-off) case. The agreement between computations and the 
experimental data is quite good up to stall. Stall is not accurately computed due to a 
lack of circumferential grid resolution on the main element [47]. Stall on a flapped 
airfoil typically occurs when the separation point on the main element suddenly 
shifts forward with a small increase in angle of attack. The relatively sparse grid 
resolution along the mid-chord region of the airfoil prevents the accurate 
computation of this separation point shift. Higher grid resolution has shown to cure 
this problem [47], but would not be practical in the three-dimensional computations. 
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The computed results show a slight decrease in lift-curve slope compared to the 
experiment. Also, an offset of 0.1 in Cl is seen at zero angle of attack. 

A similar plot for the 30° flap deflection (landing) case is shown in Figure 23. 
The trends are similar to those shown for the take-off configuration. The most 
notable difference between Figure 22 and Figure 23 is the increased offset between 
experiment and computations in the landing case. The offset between the curves is 
0.3 in Cl at zero degrees angle of attack. 

Figures 24-27 present comparisons of computed and experimental pressure 
distributions. Figure 24 shows the pressure distributions for the take-off 
configuration at 0 angle of attack. The agreement is good over the upper surface of 
the main element, but on the lower surface the agreement is poor. This was an 
unexpected result because the flow on the lower surface is largely potential and is 
usually relatively easy to predict accurately. For example, the two-dimensional 
version of the potential flow code PMARC was used to compute the pressures about 
this geometry without wind tunnel walls [48]. Figure 28 shows the resulting 
comparison. The PMARC results agree very well with the INS2D-UP results. The 
fact that these two very different methods gave such similar results provides 
confidence in the computational results. Also, INS2D-UP has proven the ability to 
accurately model multi-element airfoils in many cases [20-23]. This would suggest 
two possible sources of the disagreement: an error in the experimental data exists or 
there was some difference between the experimental and computational conditions. 
The overall agreement on the flap element is much better than on the main element. 
The most noticeable difference occurs on the upper surface at the trailing edge. The 
experiment shows a sudden pressure coefficient (Cp) increase at the most aft 
pressure tap. The computed results do not show this because no trailing edge 
bluntness was modeled. Ashby has verified that trailing edge bluntness is the 
source of this disagreement [47]. The suction peaks are captured for both elements. 



The agreement appears much better for the a =10° case shown in Figure 25. 
The only visible disagreement is in the area of the suction peak on the main element, 
INS2D-UP slightly under predicted the value of maximum suction. This small 
difference could be attributed to local compressibility in the experiment which is not 
present in the computations. 

The landing configuration cases. Figures 26-27, show results similar to those 
discussed above. However, the discrepancy between the CFD and experiment is 
visible on both the upper and lower surfaces for a=0°. As with the take-off flap 
deflection, the agreement improves at the higher angle of attack. A suction spike is 
visible at the leading edge of the main element in Figure 27. This is a local 
phenomena caused by a discontinuity in the curvature of the surface grid. Figure 29 
shows a plot of the second derivatives of the airfoil surface. The closer the curvature 
discontinuity is to the suction peak, the more noticeable the pressure spike becomes. 
A similar phenomena occurs on the flap element but is less noticeable because of its 
lower magnitude. 

5.3 FULL-SPAN FLAP RESULTS 

The full-span flap case provided a logical transition point from the two- 
dimensional airfoil work to a three-dimensional study. The goal of this phase of the 
research was to gain confidence in the ability of INS3D-UP to predict high-lift flows 
while simultaneously studying the effects of spanwise grid density. It was also 
hoped that any inherent differences between INS2D-UP and INS3D-UP would be 
apparent. These comparisons could be used to gain insight into the effects of wind 
tunnel walls and three-dimensional turbulence in two-dimensional wind tunnel 
testing. A full-span flap solution required 600 iterations to converge. For this case, 
convergence was set at the point where the divergence of velocity fell below unity. 
A slightly less stringent convergence criteria was used in three dimensions due to 
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the increase in computational requirements. Convergence took 12 hours on a Cray 
C-90 supercomputer. 

The two-dimensionality of the full-span result can be seen in Figure 30. This 
figure shows the pressure distribution for three spanwise stations. The computed 
pressures are the same at every spanwise location. The pressures did not change 
with varying spanwise grid density. Even in regions of separated flow, the 
pressures did not change noticeably, and the spanwise velocity perturbations were 
so small in magnitude that the flow remained essentially two-dimensional. This can 
be attributed to several factors. Most importantly, the grid planes were identical at 
every spanwise station. No physical wing is geometrically identical along its span, 
and these slight surface variations introduce three-dimensional disturbances into the 
flow. Turbulence itself is also inherently three-dimensional in nature. In a real, 
separated flow, the presence of turbulence is enough to generate a three-dimensional 
velocity field. Because CFD requires the use of turbulence models, some of the 
three-dimensionality of turbulence is lost. This is because the models are designed 
to approximate the eddy viscosity of turbulence, not the physical fluctuations in the 
flow. The fact that inviscid walls were modeled in the computations removed 
another possible source of three-dimensional flow disturbances. A wall with a 
boundary layer that impinges on the model has the possibility of introducing 
spanwise velocity perturbations. With an inviscid wall, no boundary layer was 
present, so two-dimensional flow was maintained at the wall. 

The pressure distributions computed with INS3D-UP are compared with the 
INS2D-UP results in Figure 31. As shown, the two codes predict nearly identical 
pressure distributions for this geometry. Figure 32, shows a good comparison 
between experimental and computed pressures. Since the full-span flap results so 
closely match the two-dimensional results, it is expected that other configurations 
would show the trends discussed in Section 5.2. This demonstrated that there were 



no fundamental disagreements between the INS2D-UP and INS3D-UP codes, giving 
some confidence that INS3D-UP could predict high-lift flows. Also, the results 
that much work, such as grid studies, could be done in two-dimensions, 
greatly reducing the computational requirements. Because of these results, it was 
felt that no additional knowledge would be gained by computing additional full- 
span configurations. 

5.4 HALF-SPAN FLAP RESULTS 

Solutions for the take-off and landing configurations at ten degrees angle of 
attack and a Reynolds number of 3.7 million were computed. In most instances, the 
effect of the flap edge was the same in both flowfields, but was exaggerated in the 
landing flap deflection case. Therefore the results presented are for the landing 
configuration unless specifically stated otherwise. Results from the take-off case are 
included where they differ significantly from those of the landing case. Both of 
these solutions were iterated until the flow divergence fell below one. A run 
required approximately 550 iterations and twenty-four Cray C-90 hours to converge. 

Figure 33 provides a general view of the flowfield about the half-span flap 
wing in the landing configuration. This figure contains particle traces colored by 
velocity magnitude. The surface particles are restricted to the surface which would 
represent oil flow patterns in an experiment. The particle traces away from the 
surface are unconstrained. The blue regions represent velocities well above free 
stream, and the red areas are velocities close to free stream. This figure illustrates 
the large scale effects of the flap edge: the flap tip vortex, large amounts of spanwise 
flow, and separation lines can all be seen. Figure 34 shows the same particle traces 
for the take-off configuration. A weaker tip vortex is present and the amount of 
spanwise flow has decreased. The separation lines have also nearly disappeared. 
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The presence of the flap edge greatly modifies the lift distribution of the wing. 
Figures 35 and 36 show plots of the spanwise lift distribution for the landing and 
take-off configurations. In both cases, an abrupt change in lift is seen at the location 
of the flap edge. This change is caused by a discontinuity in the chord and camber 
distributions as the wing section changes from the unflapped to the flapped airfoil. 
For the landing case, the total section lift coefficient, referenced to the unflapped 
airfoil chord, changes by 0.3 at the flap edge. This change in circulation (lift) causes 
vorticity to be shed into the wake in the form of a tip vortex. This vortex in turn 
induces a velocity field which further modifies the lift. The induced velocity 
appears as upwash on the unflapped wing segment and as downwash on the 
flapped half of the wing. Upwash increases the local angle of attack which increases 
the lift. This explains why the lift coefficient on the unflapped side of the wing is 
much higher than the two-dimensional value at the same angle of attack 
Conversely, the lift on the flapped half of the wing is much lower than its 
corresponding two-dimensional value. The effect of the vortex varies with the 
distance from the flap edge. If the wing span were extended infinitely in both 
directions, the lift coefficients would eventually return to the two-dimensional 
values. 

The shape of the takeoff and landing lift distributions differ slightly in the 
vicinity of the flap edge. In the takeoff case (Figure 35), the lift on the flap varies 
across the span, while almost no change in lift occurs over the span in the landing 
case (Figure 36). The velocity induced by the tip vortex changes the local angle of 
attack across the span, and the sections nearest the tip are affected the most. The 
change in the angle of attack is responsible for the lift changes for the takeoff 
configuration. In the landing case, the presence of separated flow modifies the 
shape of the lift distribution. As the local angle of attack decreases near the flap 
edge, the separation point moves aft (Figure 33). The reduction in the amount of 



separated flow causes an increase in the lift generated by the section. In this case, 
the lift gained by reducing the size of the separated region is very nearly the same 
magnitude as the decrease in lift due to the reduced angle of attack. Therefore, the 
lift on the flap remains nearly constant over the span. 

The spanwise variation in lift is caused by variations in the pressure 
distributions. Figure 37 illustrates this phenomena for two stations along the 
unflapped half-span. The increased upwash nearest the flap edge (i.e. mid-span) 
has increased the effective angle of attack compared to the section near the wall, 
decreasing the pressure coefficients on the upper surface of the wing. A slight 
increase in the lower surface pressure coefficients at mid-span can also be seen. The 
most significant difference between the two pressure distributions shown in Figure 
37 occurs at the trailing edge. The section closest to the flap edge displays a slightly 
negative pressure gradient at the trailing edge, while the section near the wall shows 
continuing pressure recovery all the way to the trailing edge. This can be explained 
by the high induced velocities at the trailing edge of the mid-span section. These 
induced velocities create a region of increased dumping velocity (the velocity at 
which the flow leaving the trailing edge must adjust to), decreasing the Cp to which 
the pressure must recover. For an airfoil without a flap, the trailing edge condition 
is such that the flow tends toward freestream as it leaves the airfoil surface. In a 
multi-element airfoil, the flow leaving the trailing edge of an upstream element often 
has the tendency to approach a velocity much higher than freestream due to the high 
velocities induced by the downstream element. The unflapped portion of the wing 
is more susceptible to this phenomena because of the exposed edge created by the 
mismatch in chord at the junction of the flapped and unflapped sections. In all 
cases, the pressure distributions for sections between the mid-span and the wall fall 
between the two curves shown in Figure 37. 
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Figure 38 shows a similar plot for the flapped main element section. In this 
case the upper surface Cp's become more negative and the lower surface Cp's 
become more positive as they move away from the mid-span toward the wall. As 
with the unflapped section, the largest differences are seen at the trailing edge. No 
upturn in the pressure distribution is seen at the trailing edge of the mid-span 
section in this case. The upper surface velocities at the flapped trailing edge must 
match the velocities on the upper surface of the unflapped portion of the wing. The 
upper surface pressures must adjust so that they are continuous across the junction 
between the flapped and unflapped sections. Also, the lower surface of the flapped 
trailing edge is perpendicular to the solid surface formed by the patch grid because 
the flapped chord is five percent shorter than the unflapped chord. This solid 
surface allows a pressure difference between the upper and lower surfaces upstream 
of the trailing edge. Away from the mid-span, the pressure distribution quickly 
begins to match the distribution at the wall. 

The greatest variation of pressure distributions occurs on the flap element. 
Figure 39 shows pressure distributions for three spanwise locations; at the flap edge, 
0.1 chord from the edge, and at the wall. At the flap edge, the pressures are 
dominated by the tip vortex which completely alters the shape of the distribution. 
As shown in the plot, the leading edge suction peak is suppressed at the flap edge. 
Away from the edge, the flow accelerates around the airfoil leading edge as it moves 
away from the stagnation point. On the edge, however, the flow is not forced to 
follow a path over the airfoil surface: it can move spanwise as well. In this case, the 
flow does move in a spanwise direction toward the edge as it leaves the stagnation 
point, relieving the suction peak. The magnitude of the relief varies with spanwise 
position. As seen in Figure 39, the peak has nearly reached the wall value at a 
distance of 0.1 chord from the edge. The velocities over the aft part of the airfoil’s 
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upper surface are higher at the edge than near the wall. The tip vortex travels just 
above the upper surface, inducing a region of high speed flow. 

A consequence of the spanwise pressure variations is a shifting of the 
separation point location along the span. In Figure 40 surface particle traces 
illustrate a separation line on the flap. The separation point is farthest forward at the 
wall and farthest aft at the tip. The separation location is closely related to the 
pressure distributions discussed in the previous paragraph. As was shown in Figure 
39, the flap pressure recovery is steepest at the wall and almost non-existent at the 
flap edge. Also shown in Figure 40 is the location of the two-dimensional separation 
line for the flapped airfoil at the same angle of attack. At the wall, the two- and 
three-dimensional flows separate very close to the same chordwise position. The 
separation points are essentially the same due to the resemblance of the 
corresponding flap pressure distributions. Figure 41 shows the half-span flap 
pressure distribution at the wall compared to the two-dimensional pressures. Since 
the distributions are almost identical, it is not surprising that the separation points 
are in the same location. The total lift generated by these sections is different (Figure 
35), which implies that the lift coefficients on the main elements must also be 
different. The different pressure distributions (different lift coefficients) would 
cause the wakes leaving the main element to be different in both cases. Since the 
flap pressures and separation points are nearly the same, it can be concluded that 
the flap boundary layer in this case is not sensitive to small changes in the main 
element wake. 

Figure 42 shows a comparison of the two and three-dimensional separation 
lines for the unflapped half-span. For this portion of the wing, the two-dimensional 
separation point is farther aft than the three-dimensional separation point at the 
wall. This is an expected result because the three-dimensional section is operating at 
a higher lift-coefficient than the two-dimensional section due to the upwash induced 
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by the tip vortex. For a given geometry higher lift usually results in earlier 
separation. Figure 43 shows a comparison of the wing surface pressures at the wall 
with the airfoil pressures. The steeper upper surface pressure recovery of the three- 
dimensional case is clearly the cause of the earlier boundary layer separation. Near 
the mid-span, however, the three-dimensional separation point is aft of the two- 
dimensional point, even though the three-dimensional section is generating more 
lift. The favorable pressure gradient at the trailing edge of the mid-span section, 
seen in Figure 37, enables the boundary layer to remain attached longer than in the 
two-dimensional case. Very close to the mid-span no clear separation line can be 
seen, as there is a large region of spanwise flow. 

Figures 44-46 show the spanwise variation of pressure distributions for the 
take-off configuration. Figures 44 and 45, which show the pressures for the main 
element segments, are very similar to Figures 37 and 38. The differences across the 
span are less apparent for the take-off configuration. This is expected because the 
trailing vortex is weaker for the take-off case than it is for the landing case: the 
vortex strength is proportional to the lift on the flap which is higher for the landing 
configuration. The flap pressures do differ between the take-off and landing 
configurations. Figure 46 shows the variation of pressure coefficients over the span 
of the flap. Not only is the shape of the distribution at the edge different, but the 
effect of the edge is felt farther away from the edge in the landing case. The 
pressures 0.1 chord from the edge in the take-off case resemble the wall pressures 
more closely than those at the wall in the landing case. This is expected because of 
the weaker tip vortex in the take-off case. The leading edge pressures are 
suppressed more in the take-off configuration. The pressure coefficients decrease 
over the aft two thirds of the flap due to the presence of the tip vortex above the 
surface. As in the landing case, the vortex induces a high local velocity near the flap 



edge, reducing the pressure. From Figure 46 it appears that most of the lift on the 
flap near the edge is vortex induced. 

The tip vortex primarily responsible for the differences between the two and 
three-dimensional flows is created by the sudden change in circulation that occurs at 
the flap edge. Physically, the vortex arises from the pressure imbalance between the 
upper and lower surfaces of the flap. High pressure flow on the lower surface of the 
flap moves around the flap edge toward the low pressure region on the upper 
surface. This rotational flow structure is converted downstream by the local flow 
velocity which turns the structure into a vortex filament. Since the pressure 
difference between the upper and lower surfaces imparted angular momentum to 
the vortex, it has a tendency to continue to rotate until it is damped out by the fluid 
viscosity. This vortex decay usually occurs far downstream of the wing tip which 
created the vortex. 

Figure 47 illustrates the location of eight positions over the upper surface of 
the flap investigated using particle traces. Figures 48-55 show the particle traces for 
these eight stations. The particles are restricted to the planes shown, so the plots 
contain projections of the streamlines onto the planes. In Figure 48 the vortex is not 
yet apparent, but the flow can be seen to travel parallel to the flap edge and around 
the comer onto the flap surface. The first signs of a vortex core are visible in Figure 
49. The vortex filament is seen to increase in size as it progresses downstream 
(Figures 50-55), and it remains above the flap from its beginning to the trailing edge. 
The flow patterns shown in Figures 48-55 correspond to the surface flow patterns 
shown in Figure 40. The flow on the surface moves away from the flap edge before 
the vortex core is present. After the vortex appears above the surface, it induces a 
velocity beneath it towards the flap edge. The particle traces in Figure 40 show the 
changing direction of the flow near the flap edge due to the vortex. The vortex for 
the take-off flap deflection travels along a similar path. 



It appears that leading edge flow near the flap edge moves in the span wise 
direction toward the flap edge. The flow separates at the flap edge forming a small 
recirculation region. The blue particle traces near the leading edge in Figure 56 
show the recirculation region. Fluid particles not caught in the recirculation region 
turn so that they flow along the flap edge. The particles begin to move towards the 
upper surface as they feel the influence of the flow leaving the lower surface. This 
can be seen in Figure 56. When they reach the upper surface, the particles form 
another small separation bubble as they flow around the comer formed by the upper 
surface and the flap edge. It is believed that this high vorticity bubble is the 
beginning of the vortex core. The bubble is barely visible at the comer in Figure 47. 

An aspect of the flow that demonstrates the resolution provided by the flap 
edge grid is the secondary recirculation regions seen in Figure 57. These flow 
structures are fed by the fluid leaving the upper and lower surfaces at the flap edge. 
The flap edge grid did a good job of resolving the flow in this region, but its small 
size made the interpretation of the results difficult. Much of the interesting flow 
physics occurs across the boundaries of the flap edge grid, requiring the tracking of 
the flow from one grid zone to another. Since the grid resolutions are somewhat 
different in each zone, smooth transitions between zones is not always possible. The 
other drawback with the current flap edge grid is that a large number of Chimera 
boundaries occur in the region of highest flow activity. The information is passed 
between zones via a low-order interpolation scheme, so accuracy is lost in regions 
where the highest accuracy is desired. This also slows convergence somewhat 
which increases the required computer time for a given case. Ideally, the flap 
surface and edge would be modeled with a single zone. 

Insufficient resolution has been identified in two regions of the grid. The 
vortex formation and dissipation were not adequately modeled by the current grid. 
The formation of the vortex is not well resolved due to a lack of spanwise and 
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chordwise grid density in the flap grid at this point. The spanwise spacing was kept 
at the current size so that adequate overlap could be obtained while keeping the 
overall grid size at a minimum. The chordwise spacing was determined from two- 
dimensional results, and the mid-chord region does not typically require many 
points to resolve for unseparated flows. In the current computations the trailing 
vortex is dissipated too quickly behind the wing. Accelerated decay occurs as the 
vortex travels into a portion of the grid that is not well resolved. The increasing cell 
size adds numerical dissipation to the solution which causes the vortex to damp out 
in a physically unrealistic fashion. Repositioning of the main element wake cut so 
that it follows the vortex path would better resolve the flow and may maintain the 
vortex in a more realistic manner. These grid modifications were not performed as 
part of this study. 



CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 

A numerical investigation of a simple high-lift system was performed using 
an incompressible Navier-Stokes flow solver. The study began by investigating 
take-off and landing configurations of a two-element, two-dimensional airfoil at a 
range of angles of attack. Grid issues were addressed in two dimensions due to the 
much lower computational requirements. The two-dimensional grid planes were 
converted to a three-dimensional, full-span flap grid. This served as a validation 
case for the three-dimensional flow solver. Once confidence in the three- 
dimensional flow solver was established, a grid scheme was developed to model the 
effects of a half-span flap. The half-span flap case provided a means to study the 
flow physics associated with the flap edge. 

The full-span results agreed well with experimental data for most cases. 
Some differences in the pressure distributions were seen between the computed and 
experimental results, particularly on the lower surface of the airfoil with high flap 
deflections. Because three computational techniques (INS2D-UP, INS3D-UP, and 
PMARC) produced very similar results, it is believed that the disagreement lies in an 
experimental error or in a modeling inconsistency between the experiment and the 
computations. 

The grid scheme developed to represent the half-span flap configuration was 
successful. The grid provided enough resolution to capture the large-scale influence 
of the trailing vortex on the lift distribution of the wing. The change in the pressure 
distributions created large amounts on spanwise flow and modified the separation 
lines compared to the two-dimensional case. Secondary recirculation regions were 
seen in the solution at the comers formed between the flap surface and the flap edge. 
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The details of the tip vortex formation were not clearly seen in the current results, 
but the vortex trajectory could be followed over the upper surface of the flap once it 
has formed. 

Findings from this research provided a number of possible extensions of this 
study. The effects of grid density at the flap edge could be investigated. An increase 
in the number of spanwise and chordwise points at the edge may show the 
formation of the tip vortex in more detail. Also, creating a single zone grid around 
the flap edge would better resolve the flow in this area without any uncertainty due 
to the Chimera boundaries. The wake cut of the main element could also be 
positioned so that it follows the trajectory of the tip vortex. This would make use of 
the large number of points along the wake cut and would more realistically model 
the convection of the trailing vortex. 

Steady state computations were performed in this study to reduce the 
computational requirements of the research. Near the flap edge, a number of 
unsteady phenomena, such as separated flows, exist. Therefore a time-accurate 
investigation may be required to truly model the flow in this region. An unsteady 
analysis is certainly required to track any movement of the vortex core, which could 
influence the rest of the flowfield. 
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Figure 1: Trends in Boeing transport high-lift system development 

(from Reference 3). 
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Figure 2: Lift efficiency versus high-lift system complexity (from Reference 4) 





Figure 3: Curved computational mesh. 



Figure 4: Relationship between physical and computational domains. 



Figure 5: NACA 632-215 Mod. B airfoil. 



Figure 6: NACA 632-215 Mod. B airfoil with 30% chord Fowler flap 

in take-off configuration. 
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Figure 7: NACA 632-215 Mod. B airfoil with 30% chord Fowler flap 



Figure 8: Definition of gap and overlap. 
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Figure 10: A two-dimensional grid plane. 



Figure 11: Overset grids before removing points inside flap surface. 
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Figure 13: Use of hole-cuts to create overset grids. 
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Figure 16: Landing configuration airfoil 


Patch grid 


Figure 17: Exposed surface at flapped-unflapped intersection. 



Figure 18: Patch grid plane. 









Angla of Attack 

Figure 22: Lift coefficient versus angle of attack for two-dimensional take-off 

configuration. 



Angle of Attack 

Figure 23: Lift coefficient versus angle of attack for two-dimensional landing 

configuration. 
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Figure 24: Pressure distribution for two-dimensional take-off configuration, a=0°. 
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Figure 25: Pressure distribution for two-dimensional take-off configuration, a =10°. 
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Figure 27: Pressure distribution for two-dimensional landing configuration, a =10 
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Figure 30: Pressure distributions for three spanwise stations on full-span flap, 

a= 10°. 
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Figure 31: Pressure distributions computed with INS2D-UP and INS3D-UP for full- 

span flap, a- 10°. 










Figure 33: Particle traces colored by velocity magnitude for half-span flap in landing 

configuration, a =10°. 
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Figure 34: Particle traces colored by velocity magnitude for half-span flap in take-off 

configuration, a- 10°. 
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Figure 37: Spanwise pressure variation for landing configuration unflapped main 

element, a =10°. 
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Figure 38: Spanwise pressure variation for landing configuration flapped main 

element, a =10°. 
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Figure 41: Pressure distribution comparison between half-span and airfoil 
computations for landing configuration flap element, a -10°. 
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Figure 43: Pressure distribution comparison between half-span and airfoil 
computations for landing configuration unflapped element, a =10°. 
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Figure 44: Spanwise pressure variation for landing configuration unflapped main 

element, a =10°. 
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Figure 45: Spanwise pressure variation for landing configuration flapped main 

element, a =10°. 



Figure 46: Spanwise pressure variation for landing configuration flap element, 

a =10°. 
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Figure 55: Particle traces at station 8 on flap element in landing position. 



Figure 56: Particle traces colored by pressure at the flap edge of landing 

configuration. 
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